Space-time wave packets localized in all dimensions

Optical wave packets that are localized in space and time, but nevertheless overcome diffraction and travel rigidly in free space, are a long sought-after field structure with applications ranging from microscopy and remote sensing, to nonlinear and quantum optics. However, synthesizing such wave packets requires introducing non-differentiable angular dispersion with high spectral precision in two transverse dimensions, a capability that has eluded optics to date. Here, we describe an experimental strategy capable of sculpting the spatio-temporal spectrum of a generic pulsed beam by introducing arbitrary radial chirp via two-dimensional conformal coordinate transformations of the spectrally resolved field. This procedure yields propagation-invariant ‘space-time’ wave packets localized in all dimensions, with tunable group velocity in the range from 0.7c to 1.8c in free space, and endowed with prescribed orbital angular momentum. By providing unprecedented flexibility in sculpting the three-dimensional structure of pulsed optical fields, our experimental strategy promises to be a versatile platform for the emerging enterprise of space-time optics.


FACE OF THE LIGHT-CONE
In our previous work on space-time (ST) wave packets in the form of light-sheets [1][2][3][4][5][6], we make heavy use of the representation of the wave-packet spectral support domain on the surface of the light-cone. This is a useful visualization tool that provides physical intuition with regards to the structure and behavior of ST wave packets. In this Section, we briefly review this representation in the reduced-dimension case of ST light sheets [2], where the field is localized along one transverse dimension and extended uniformly along the other. We refer to these field structures as 2D ST wave packets (one transverse dimension and one longitudinal dimension).
We then proceed to show that such a representation can also be gainfully employed with minor changes for ST wave packets localized in all dimensions. We refer to such field structures as 3D ST wave packets (two transverse dimensions and one longitudinal dimensions). Therefore, the wealth of results that have amassed over the past few years based on this conceptual framework [6] can be appropriated for the new 3D ST wave packets investigated here.

A. Light-cone representation for 2D ST wave packets (light sheets)
When the field is held uniform along one transverse dimension (say y), then the dispersion relationship in free space is k 2 x +k 2 z =( ω c ) 2 , where k x and k z are the transverse and longitudinal components of the wave vector along x and z, respectively, ω is the temporal frequency, and c is the speed of light in vacuum. This relationship is represented geometrically by the surface of a cone that we refer to as the light-cone. A monochromatic plane wave e i(k x x+k z z−ωt) is represented by a point on the surface of the light-cone ( Supplementary Fig. 1). In general, a pulsed beam E(x,z;t) is expressed as a product of a slowly varying envelope ψ(x,z;t) and a carrier term where ω o is a fixed temporal frequency, and k o =ω o /c is its associated wave number.
The envelope is written in terms of an angular spectrum as follows: ψ(x,z;t)= dk x dΩ ψ(k x ,Ω)e i{k x x+(k z −k o )z−Ωt} , (S1) where Ω=ω −ω o , and the spatio-temporal spectrum ψ(k x ,Ω) is the 2D Fourier transform of ψ(x,0;t). The spectral support domain for a pulsed beam or wave packet corresponds in general to a 2D area on the surface of the light-cone [2]; see Supplementary x +k 2 z =( ω c ) 2 with a tilted spectral plane in the region k z >0. The spectral projection onto the (k z , ω c )-plane is a straight line, and that onto the (k x , ω c )-plane is a segment of a conic section. (b) The spatio-temporal intensity profile in the plane z=0, I(x,z=0;t), and the transverse intensity profile in the same z=0 plane at the wave-packet center t=0, and off-center t>0. In the (x,t)-plane, the spatio-temporal intensity profile is X-shaped. In the (x,y)-plane, the 2D ST wave packet at t=0 takes the form of a light sheet that is localized along the x-axis and extends uniformly along the y-axis.
For the special case of a propagation-invariant 2D ST wave packet, the spectral support domain is confined to the intersection of the light-cone with a spectral plane that is parallel to the k x -axis and makes an angle θ with the k z -axis, which is thus given by the equation: where we refer to θ as the spectral tilt angle. Consequently, there is a one-to-one relationship between Ω and |k x |, which takes the form of a parabola in the narrowband (∆ω ω o , where ∆ω is the temporal bandwidth) and paraxial (∆k x k o , where ∆k x is the spatial bandwidth) limits: where n=cotθ is the wave packet group index in free space. Because of this correspondence between Ω and |k x |, the spatio-temporal spectrum has a reduced dimensionality with respect to a conventional pulsed beam ψ(k x ,Ω)→ ψ(k x )δ(Ω−Ω(k x )), where Ω=Ω(k x ) is given by Eq. S3.
The wave-packet envelope now takes the simpler integral form: which is a 2D ST wave packet that travels rigidly in free space at a group velocity v=ctanθ.
The spectral projection of the 2D ST wave packet onto the (k z , ω c )-plane is a straight line that makes an angle θ with the k z -axis and intersects with the light-line k z = ω c at the point (k z , ω c )=(k o ,k o ). The corresponding spectral projection onto the (k x , ω c )-plane is a segment of a conic section: an ellipse when 0 • <θ<45 • or 135 • <θ<180 • , whereupon | v|<c; a hyperbola when 45 • <θ<135 • , whereupon | v|>c; a straight tangent line when θ=45 • and v=c, corresponding to a plane-wave pulse; and a parabola when θ=135 • and v=−c. We refer to ST wave packets associated with the range 0 • <θ<45 • as subluminal, with 45 • <θ<90 • as superluminal, and with θ>90 • (whereupon v<0) as negative-v ST wave packets.
Note that causal emission and propagation require that only the values k z >0 be considered, so that the light-cone half corresponding to the acausal backward-propagating components k z <0 is eliminated from consideration [4,7].
The one-to-one correspondence between |k x | and ω has a crucial impact on the form of the axial evolution of the time-averaged intensity I(x,z)= dt|E(x,z;t)| 2 = dt|ψ(x,z;t)| 2 . Substituting for ψ(x,z;t) from Eq. S4 we reach: where . In other words, I(x,z) is altogether independent of the axial coordinate z, and is formed of the sum of a constant background pedestal term I o and a localized spatial feature at x=0. Moreover, the height of the localized spatial feature cannot exceed the height of the pedestal. This structure of the intensity profile for 2D ST wave packets has been borne out in previous measurements [2,5].
It is important to note that this structure is unique to 2D ST wave packets. The time-averaged intensity of 3D ST wave packets cannot be separated into a sum of a pedestal and a localized central feature (see main text).
B. Spectral representation for 3D fields on the light-cone surface

Conventional optical fields
When both transverse coordinates x and y are retained, we have the general dispersion relationship k 2 x +k 2 y +k 2 z =( ω c ) 2 in free space. This relationship is represented mathematically by a hypercone in 4D, which cannot be visualized in 3D space. However, because we are mostly interested in cylindrically symmetric fields, we can write the dispersion relationship as x +k 2 y is the radial wave number. This dispersion relationship can indeed be represented in 3D space. Because k r is positive-valued only, in contrast to k x in the case of the 2D ST wave packets that can take on either positive or negative values, only the quarter of the light-cone corresponding to k r >0 and k z >0 need be retained here.
A wave packet is once again given in terms of a carrier term and a slowly varying envelope E(x,y,z;t)=Re e i(k o z−ω o t) ψ(x,y,z;t) , where: and ψ(k x ,k y ,Ω) is the 3D Fourier transform of ψ(x,y,0;t). We switch to transverse polar coordinates (r, ϕ) in physical space and to the corresponding polar coordinates (k r ,χ) in Fourier space.
In physical space we have the relationships: and in Fourier space we have We can thus rewrite the angular spectrum of the envelope as follows: ψ(r, ϕ,z;t)= dk r dχdΩ k r ψ(k r ,χ,Ω) e ik r rcos(ϕ−χ) e i(k z −k o )z e −iΩt ; where the integral over χ extends from 0 to 2π, that over k r extends from 0 to ∞, and that for Ω from −∆ω/2 to ∆ω/2.
We can separate the spatio-temporal spectrum ψ(k r ,χ,Ω) with respect to the radial and azimuthal coordinates k r and ξ, respectively, as follows: whereupon the wave packet envelope can be expressed as: where we have made use of the identity 2π J (x)= π −π dye i(xsiny− y) , and J (·) is the th -order Bessel function of the first kind.
One simplification is obtained by assuming that the spatio-temporal spectrum is azimuthally symmetric, which corresponds to setting =0, and the summation over is removed. This assumption is equivalent to setting ψ(k r ,χ,Ω)→ ψ(k r ,Ω), which is thus independent of χ, whereupon: ψ(r, ϕ,z;t)= dk r dΩ k r ψ(k r ,Ω)J 0 (k r r) e i(k z −k o )z e −iΩt =ψ(r,z;t), where J 0 (·) is the zeroth-order Bessel function of the first kind. As a result of the azimuthal symmetry of the spectrum, the wave-packet envelope is also azimuthally symmetric in physical space. We will address shortly the more general case of fields with azimuthal variation.
Comparing Eq. S1 for ψ(x,z;t) to Eq. S12 for ψ(r,z;t), we find that both have 2D spatio-temporal spectra: ψ(k x ,Ω) for the former and ψ(k r ,Ω) for the latter. The light-cone in (k r ,k z , ω c )-space can thus be used to represent the spectral support domain for ψ(r,z;t). Here, each point on the light-cone at coordinates (k r ,k z , ω c ) corresponds to a monochromatic Bessel beam E(r,z;t)= J 0 (k r r)e i(k z z−ωt) rather than a monochromatic plane wave; see Supplementary Fig. 2.
In the case of light sheets in which the field is uniform along y, the light-cone in Supplementary   Fig. 1(a) suffices to capture the complete description of the spectral support domain. On the other hand, in the case of fields in 3D space, the light-cone in Supplementary Fig. 2(a) does not convey the whole picture because we collapsed all the plane waves having spatial-frequency pairs (k x ,k y ) into their radial counterpart k r . The picture can be completed by adding a second spectral representation in (k x ,k y , ω c )-space. Each point with coordinates (k x ,k z , ω c ) corresponds to a monochromatic plane wave e i(k x x+k y y+k z z−ωt) , where k z = ( ω c ) 2 −k 2 x −k 2 y . Although k z is not Supplementary Fig. 2. Representation of the spectral support domain of a monochromatic Bessel beam in 3D space on the surface of the light-cone. (a) A point on the surface of the light-cone k 2 r +k 2 z =( ω c ) 2 at coordinates (k r ,k z , ω c ) corresponds to a monochromatic Bessel beam of the form E(r,z;t)=J 0 (k r r)e i(k z z−ωt) . (b) In (k x ,k y , ω c )-space, the point on the light-cone in (a) takes the form of a horizontal iso-frequency circle of radius k r at a height ω c . (c) In physical space at an axial plane z=0, the intensity I(x,y,0;t) is uniform along t and takes the form of a Bessel function J 2 0 (k r r) in the (x,y)-plane at any t.
represented in this space explicitly, it can nevertheless be found by referring to the light-cone in Supplementary Fig. 2(a). The single point representing a monochromatic Bessel beam on the lightcone in Supplementary Fig. 2(a) corresponds to a horizontal circle of radius k r in (k x ,k y , ω c )-space in Supplementary Fig. 2 The representation in (k x ,k y , ω c )-space is particularly useful for wave packets in 3D space because it provides the spatio-temporal spectral structure required to synthesize the field in question. For a monochromatic Bessel beam, Supplementary Fig. 2(b) points to the well-known approach for producing such a beam by inserting a spatial filter in the Fourier domain in the form of a thin annulus, followed by a spherical converging lens [8]. (b) In (k x ,k y , ω c )-space, the spectral support domain is a horizontal iso-frequency disc at ω=ω o . The disc radius corresponds to the maximum spatial frequency in (a). (c) In physical space at z=0, the intensity I(x,y,0;t) is uniform along t with the transverse beam profile given by the plot in the (x,y)-plane at any t.
More generally, the spectral support domain for a monochromatic beam at frequency ω o is restricted to the circle at the intersection of the light-cone with a horizontal iso-frequency plane  transverse features, and whose longitudinal and temporal shape is determined exclusively by ψ(Ω), therefore: The spectral support domain in (k x ,k y , ω c )-space lies along the ω c -axis, with k x =k y =0 [Supplementary Fig. 4 Finally, the most general conventional azimuthally symmetric pulsed beam or wave packet has a 2D spectral support domain on the surface of the light-cone in (k r ,k z , ω c )-space as shown in Supplementary Fig. 5(a). Typically, ψ(k r ,Ω) is separable with respect to k r and Ω, ψ(k r ,Ω)→ takes the form approximately of a cylindrical volume centered on the ω c -axis, whose radius is the maximum width of spatial spectrum ψ r (k r ), and whose heights is the extent of the temporal spectrum ψ t (Ω).

3D ST wave packets
When considering 3D ST wave packets that are localized along all dimensions, such that the time-averaged intensity takes the form of an axial needle [9][10][11], the spectral support domain is once again restricted to the intersection of the light-cone k 2 r +k 2 z =( ω c ) 2 with the spectral plane Ω=(k z −k o )ctanθ, as in the case of the 2D ST wave packets, where the plane is parallel to the Supplementary Fig. 6. Representation of the spectral support domain for a superluminal 3D ST wave packet on the surface of the light-cone. (a) The spectral support domain is restricted to the hyperbola on the light-cone surface at its intersection with a spectral plane that is parallel to the k r -axis and makes an angle θ with the k z -axis (45 • <θ<90 • ). (b) In (k x ,k y , ω c )-space, the spectral support domain is a 2D surface in the form of one half of a two-sheet hyperboloid (an elliptic hyperboloid). (c) In physical space at z=0, the intensity I(x,y,0;t) takes the form of two cones emanating from the central feature at t=0. In any meridional plane such as x=0, the spatio-temporal intensity profile is X-shaped. In the (x,y)-plane at t=0, the intensity profile is centered at x=y=0. At t =0, the profile in the (x,y)-plane is an annulus whose radius increases with t. k r -axis and makes an angle θ with the k z -axis. The group velocity of the propagation-invariant 3D ST wave packet is v=ctanθ: The impact of θ on the shape of the conic section resulting from this intersection is the same as that for the 2D ST wave packet or light sheets. We depict the superluminal scenario in Supplementary has an ellipse for its spectral support domain. In the narrowband paraxial regime, the conic section can be approximated in the vicinity of k r =0 by a parabola: . (S15) Introducing this quadratic relationship between Ω and k r into a pulsed beam is the goal of our synthesis methodology. Note that for small bandwidths, this also entails a quadratic relationship between λ−λ o and k r , where λ o = 2π k o . In the case of the 3D ST wave packets, the structure of spectral representation in (k x ,k y , ω c )space is particularly instructive, as shown in Supplementary Fig. 6 Supplementary Fig. 6(b). Because each k r is associated with a different ω, the circles in Supplementary Fig. 6(b) are all located at different heights, thus forming a 2D surface rather than the 3D volume in Supplementary Fig. 5(b) for a conventional pulsed beam. Because the spectral trajectory on the light-cone surface in Supplementary Fig. 6(a) is a hyperbola, the surface in Supplementary Fig. 6(b) is one half of a two-sheet hyperboloid (an elliptic hyperboloid) centered on the ω c -axis. In physical space, as shown in Supplementary Fig. 6(c), the spatio-temporal intensity profile at a fixed axial profile z=0 takes the form of a central peak at x=y=0 and t=0, with two cones centered on the t-axis emanating from this peak. Consequently, in any meridional plane, for example x=0, the intensity profile is X-shaped, similarly to its 2D ST wave packet counterpart [ Supplementary Fig. 1(c)]. The shape of the cross section of the intensity profile is time-dependent: at t=0 it is localized at x=y=0, at t =0 it takes the form of an annulus whose radius increases with t.
The corresponding graphs for a subluminal 3D ST wave packet are plotted in Supplementary   Fig. 7(b,c). Because the intersection of the spectral plane with the light-cone is an ellipse, the spectral support domain in (k x ,k y , ω c )-space is an ellipsoid of revolution, or spheroid, as shown in Supplementary Fig. 7(b). Depending on the value of θ, this spheroid may be oblate of prolate. of the spectral support domain on the light-cone surface onto the (k r , ω c )-plane is a circle, and the spectral support domain in (k x ,k y , ω c )-space is a sphere. Of course, in all cases only the plane-wave components corresponding to k z >0 are physically meaningful. The dynamics of 3D ST wave packets upon free space propagation are calculated by applying Fresnel propagation to the wave packets after introducing the space-time coupling given in Eq. S15 into the optical field.

Supplementary Note 2. SYNTHESIS OF 3D SPACE-TIME WAVE PACKETS
The experimental methodology to synthesize 3D ST wave packets shown in Fig. 3 in the main text is expanded in more technical detail here in Supplementary Fig. 8. Our strategy consists of three stages as outlined in Supplementary Fig. 8 We start off with femtosecond plane-wave pulses from a mode-locked Ti:sapphire laser (Tsunami; Spectra Physics) of width ≈100 fs and bandwidth ∆λ≈10 nm centered at a wavelength of ≈800 nm. The pulses are directed to the first stage of the synthesis system: spectral analysis using a volume chirped Bragg grating (CBG), which spatially resolves the spectrum [12,13]. A double-pass through the CBG produces a linear spatial chirp but with a flat phase front.
The spectrally resolved wave front from the CBG arrangement is then fed to the second stage: a 1D conformal mapping implemented by a pair of spatial light modulators (SLM; Meadowlark 1920×1080 series) that produces a logarithmic coordinate transformation to yield a logarithmic spatial chirp. The transformed wave front is then directed to a fixed log-polar-to-Cartesian coordinate transformation. This 2D coordinate transformation maps a line at its input into a circle at its output [14][15][16], and is implemented by means of two refractive [17] or diffractive [18] phase plates. The combination of the 1D spectral transformation and the 2D coordinate transformation produces a field endowed with a quadratic radial chirp. Finally, a spherical converging lens performs an optical Fourier transform along both transverse dimensions to yield 3D ST wave packets. We proceed to provide a detailed description of each stage of this novel spatio-temporal synthesis setup.

Introducing spatial chirp
The goal of the first stage in the setup is to spatially resolve the spectrum but retain a flat phasefront, which is necessary for the successful operation of the subsequent coordinate transformations. spectrum does not have a flat phase. Instead, we make use of an arrangement based on a volume CBG to achieve this goal.
The CBG is a reflective Bragg grating [19] with a multilayered structure having a linearly varying periodicity Λ(z) along the longitudinal axis z. Consequently, different wavelengths are reflected from different depths z within the grating volume [13]. As a result, the CBG introduces a spectral chirp into normally incident pulses, thus stretching the plane-wave pulse in time [ Supplementary Fig. 9(a)]. For this reason, volume CBGs are widely used in high-power chirped pulse amplification (CPA) systems, where they are well-known for their high damage threshold and their ability to introduce extremely large spectral chirp [20,21].
Our goal is to produce spatial chirp rather than spectral chirp. At oblique incidence, however, both spectral and spatial chirps are introduced; i.e., the pulse is stretched longitudinally in time, and the spectrum is resolved transversely in space [ Supplementary Fig. 9(b)]. The spatial chirp can be determined easily from two parameters: (1) the chirp rate of the CBG structure β= dλ dz = 1 2n dΛ dz , where n is the average refractive index of the CBG; and (2) the incident angle φ 1 with respect to the normal to the grating structure. After reflecting from the CBG at oblique incidence, the spectrum is spatially resolved along the x-axis such that each wavelength λ is located at a position x(λ) given by: where ζ(φ 1 )= nsin2φ 1 For our purposes here, we aim at retaining the spatial chirp while eliminating the accompanying spectral chirp [12], which we achieve by directing the field to an identical CBG placed in a reversed geometry with respect to the first one [13]; see Supplementary Fig. 9(c). The output from CBG 1 first passes through a 4 f imaging system that flips the field along x and thus reverses the sign of the spatial chirp. Consequently, after passing through CBG 2 with an opposite sign of chirp β 2 =−β 1 , CBG 2 doubles the the spatial chirp, while cancelling out the spectral chirp introduced by CBG 1 . As a result, we obtain a spatially resolved spectrum with a flat phase front.
In our setup, we used a folded configuration in which the beam is first incident obliquely on one port of the CBG, and the reflected and flipped field is then directed to the second port of the size upon incidence on the second port. This double-pass configuration produces a spatially resolved spectrum with the following distribution: where α= 1 β ζ(φ 1 ) is the linear spatial chirp rate.
Supplementary Fig. 10. The measured spatial spread of the spectrum after the spectral analysis stage in Supplementary Fig. 8(b). The measurements confirm a linear spatial chirp over a temporal bandwidth of ∆λ≈0.8 nm spread over ≈18 mm in space horizontally along x 1 . The straight line is a theoretical fit, and symbols are data points. Error bars correspond to the spectral resolution of the Optical Spectrum Analyzer (OSA) used in the measurements, δλ=10 pm.
In our experiments, we made use of a CBG (OptiGrate L1-021) with a central periodicity of Λ o =270 nm, a chirp rate of β=−30 pm/mm, an average refractive index of n=1.5, and a length of L=34 mm. The input beam is incident at an angle φ 1 ≈16 • with respect to the normal to the CBG entrance surface. We measure the spatially resolved spectrum by scanning a single-mode fiber (Thorlabs 780HP) connected to an optical spectrum analyzer (OSA; Advantest AQ6317B).
The measured spectrum after the double-pass configuration is plotted in Supplementary Fig. 10, which verifies that the CBG produces a linear spatial chirp of rate α=−22.2 mm/nm (from Eq. S17) centered at λ o =796.1 nm. Therefore, at the output of the spectral-analysis stage we have a collimated, spectrally resolved optical field.

Spectral uncertainty induced by CBGs
A critical parameter characterizing ST wave packets in general is the so-called 'spectral uncertainty' δλ: the unavoidable 'fuzziness' in the association between spatial frequencies (k x for 2D ST wave packets and k r for 3D ST wave packets) and the wavelength λ. In the case of a CBG, its spectral resolving power (the smallest spectral shift that can be resolved) determines δλ [22].
In the case of surface gratings, the spectral uncertainty at a central wavelength λ o is given by δλ=λ o /N, where N is the number of grating rulings covered by the incident beam. Therefore, δλ here depends on the beam size and the grating period. To minimize δλ of a surface grating with a fixed ruling density, it is therefore desirable to use the largest possible incident beam size. In contrast, the spectral resolving power of a CBG at normal incidence δλ BG -and thus the associated spectral uncertainty -is determined by its refractive-index modulation contrast and the grating length [13]. Furthermore, an additional contribution to the spectral uncertainty δλ GE arises at oblique incidence due to a geometric effect stemming from the finite spatial width of the input beam. When a beam of finite transverse width w is obliquely incident on the CBG, the shifted spectral intervals reflected from different points within the CBG will overlap, thereby leading to an additional contribution to the spectral uncertainty, which is proportional to the input beam width, δλ GE ∝w. The total spectral uncertainty of the CBG at oblique incidence is estimated to be δλ CBG = (δλ BG ) 2 +(δλ GE ) 2 . When the beam passes through the CBG twice, a factor-of-2 drop is expected in the spectral uncertainty due to the doubling of the total CBG length traversed by the beam.
In our experiments, we measured δλ CBG as a function of the input beam width w after traversing the CBG once [ Supplementary Fig. 11(a)] and twice [ Supplementary Fig. 11(b)]. The beam width w is tuned using a variable-width aperture preceding the CBG [ Supplementary Fig. 8(b)].
The measurements confirm that δλ CBG increases with w as expected. However, this trend ends when w∼2 mm. Below this beam size, the spectral uncertainty begins to increase rather than to decrease further, thereby leading to a valley in the plots in Supplementary Fig. 11(a,b) in the vicinity of w∼2 mm. This unexpected increase in δ CBG at small w is most likely caused by diffraction inside the CBG resulting from the small input beam size. The optimum beam width is thus w≈2 mm, which yields a minimum spectral uncertainty of δλ CBG ≈35 pm in the double-pass configuration. The synthesis experiments we performed made use of w≈2 mm, and thus we take δλ≥35 pm for the 3D ST wave packets produced.

B. Tunable 1D spectral transformation
The second stage of the 3D ST wave-packet synthesis setup is the tunable 1D spectral transformation. This is a coordinate transformation performed along the horizontal x-axis. Because the wavelengths are arranged linearly along x after the CBG and the field is uniform along y, Supplementary Fig. 11. Spectral uncertainty for the CBG. The x-axis at the input plane is labelled x 1 and that at the output plane x 2 . The targeted transformation then takes the form: where A and B are the transformation parameters, the physical significance of which will be discussed shortly. The uniform field distribution along y remains intact. This conformal mapping is implemented by means of two phase patterns placed at the input and output planes and separated by a distance d 1 . The first phase distribution Φ 1 (x 1 ,y 1 ) at the input plane performs the particular transformation given in Eq. S18. In other words, the wavelength λ located at position x 1 is now located at position x 2 . However, such a transformation does not produce a collimated field. The second phase distribution Φ 2 (x 2 ,y 2 ) placed at the output plane collimates the transformed wave front to yield an afocal transformation. Usually, a lens is placed between the two phase plates in a 2 f -configuration. In our experiments, we do not make use of a lens and instead distribute the quadratic phase associated with such a lens between the two phase plates.
The required phase profiles Φ 1 (x 1 ,y 1 ) and Φ 2 (x 2 ,y 2 ) -both of which are independent of y -can be derived using the methodology outlined in [15], and they take the form: where k= 2π λ is the wave number. In the setup, Φ 1 (x 1 ,y 1 ) and Φ 2 (x 2 ,y 2 ) are displayed on two phase-only SLMs with d 1 =400 mm [ Supplementary Fig. 8(b)]. Because the phase profiles according to Eq.S19 depend only on x and not y, 1D SLMs can be used here in principle. However, as we show later, utilizing 2D SLMs provides the possibility of also modulating the field along y, which translates into an azimuthal modulation after the subsequent 2D coordinate transformation.
In our experiments, we fix A=0.5 mm and tune B over the range B=[−15,20] mm to control the group velocity of the 3D ST wave packets, as explained in detail later. See Supplementary   Fig. 12(a) for a plot of the relationship between x 1 and x 2 when making use of these values.
The field after this 1D spectral transformation is imaged by a one-to-one 4 f system comprising two spherical lenses L s3 and L s4 , as shown in Supplementary Fig. 8(b), from the output plane of the 1D spectral transformation (x 2 ,y 2 ) to the input plane (x 3 ,y 3 ) of the 2D coordinate transformation.
Thus, the beam is flipped along the x and y axes after this 4 f system: x 3 =−x 2 and y 3 =−y 2 . This is important to keep in mind for determining the required transformation parameters A and B based on the desired group velocity v. In addition, a spatial filter is placed in the Fourier plane of the 4 f system to eliminate the undesired zeroth-order field component resulting from the limited efficiency of SLM 1 and SLM 2 .

C. Fixed 2D coordinate transformation
The third stage of the spatio-temporal synthesis strategy for 3D ST wave packets is a fixed 2D coordinate transformation; see Supplementary Fig. 8(b) and Supplementary Fig. 13. This transformation performs a conformal mapping from log-polar to Cartesian coordinate systems.
The input plane is spanned by Cartesian coordinates (x 3 ,y 3 ) and the output by (x 4 ,y 4 ). The x 4 ); see Supplementary Fig. 12(b). The transformation is given explicitly in the following form [14,15,23]: The transformation parameter D is chosen to map the vertical size of the input field y 3 = [−y max 3 ,y max 3 ] to the range ϕ=[−π,π], and thus we impose D= y max 3 π ; the value of C is selected based on the aperture size of the optics used. The 2D coordinate transformation therefore maps a vertical line located at x 3 at the input plane into a circle of radius r at the output plane. Shifting the location of the line horizontally along x 3 at the input thus results in a radial expansion or shrinkage of the circle radius at the output according to the direction of the shift at the input [ Supplementary Fig. 13(a,b)]. If the input beam is a rectangle of width ∆x 3 , the transformed beam is an annulus of radial thickness ∆r, where the inner and outer radii of the annulus correspond to the two vertical boundaries of the input rectangle. This 2D coordinate transformation is implemented by two 2D phase patterns (both depending on the x and y coordinates) separated by a distance d 2 : Φ 3 at the input plane (x 3 ,y 3 ) and Φ 4 at the output plane (x 4 ,y 4 ). For the mapping given in Eq. S20, the phase patterns take the following form [15][16][17]: The phase plates are designed such that the aperture size is 2y max 3 =8 mm, the distance separating them in the setup is d 2 =310 mm, and the transformation parameters are C=4.77 mm, D=1 mm.
where atan2(y 4 ,x 4 ) is 2-argument arctan function. Usually, a spherical lens is placed midway between the two phase patterns in a 2 f configuration. Instead, the appropriate quadratic phases are added in Φ 3 and Φ 4 , which are the last terms in Eq. S21(a) and Eq. S21(b) [17].
As discussed in the Methods Section in the main text, the 2D transformation was implemented using diffractive optics [16,[24][25][26] and separately refractive optics [17] to imprint the phase profiles in Eq. S21 via diamond-edged refractive phase plates [17] and analog diffractive phase plates [26]. The refractive optical elements used in our experiments are similar to those outlined by Lavery et al. in [17], in which the transformation parameters are C=4.77 mm, D= 3.2 π ≈1 mm, and d 2 =310 mm. Each phase plate is made of the polymer PMMA (Poly methyl methacrylate) with accurately manufactured height profiles Z 1 (x 3 ,y 3 ) and Z 2 (x 4 ,y 4 ) to imprint the required phase profiles given in Eq. S21. The phase encountered by light at a wavelength λ traversing a height Z of a material of refractive index n -with respect to the phase encountered over the same distance in vacuum -is given by Φ=2π(n−1)Z/λ. Thus, the height profile of the first element Supplementary Fig. 14(a)] and that of the profile of the second element is Z 2 (x 4 ,y 4 )= λ 2π(n−1) Φ 4 (x 4 ,y 4 ) [ Supplementary Fig. 14(b)]. See Methods in the main text for details of the fabrication and characteristics of the fabricated elements [27]. The diffractive phase plates were fabricated in fused silica using Clemson University facilities using the process outlined in [18]. See Methods in the main text for the design parameters of the phase plates.

D. Spatial Fourier transform
After the 2D coordinate transformation, the spatially resolved wavelengths are arranged radially r(λ). A spherical converging lens L s5 [ Supplementary Fig. 8(b)] performs a 2D Fourier transform to produce the 3D ST wave packet in physical space. In effect, each spatial position is now mapped into a spatial frequency k r (λ)=k r f , where f is the focal length of the lens. because of the small bandwidth ∆λ utilized in our experiments, we can use the approximation k≈k o in the scaling of the radial spatial frequency.

E. Full system analysis
We can now combine all the stages in our synthesis system to identify the role played by the various parameters involved: the chirp rate α from the CBG; the parameters A and B from the tunable 1D spectral transformation; the parameters C and D from the fixed 2D coordinate transformation; and the focal length f of the Fourier-transforming lens. After the CBG we have x 1 (λ)=α(λ−λ o ); after the tunable 1D spectral transformation we have x 2 (λ)=Aln{ x 1 (λ) B }; after the 4 f imaging system we have x 3 =−x 2 ; after the fixed 2D coordinate transformation we have r(λ)=Cexp{− x 3 (λ) D }; and, finally, after the Fourier-transforming lens we have k r =k o r(λ) f . Combining all these steps, we have: To obtain the desired radial chirp from Eq. S15, we must impose the constraint D A =2, whereupon: The parameter C=4.77 mm is determined by the size of the optics (the aperture size is 2y max 3 = 8 mm), f =300 mm, α=−22.2 mm/nm is determined by the CBG, and λ o ≈798 nm. This leaves B as the sole free parameter determining the chirp rate, and thus also the group velocity v.
By comparing the result in Eq. S23 to the target spatio-temporal spectrum for 3D ST wave packets given in Eq. S15, we can obtain the group index: As a result, we can generate ST wave packets with group velocities from subluminal (B<0) to superluminal (B>C 2 |α|λ o /2 f 2 ), and in principle negative group velocities (0<B<C 2 |α|λ o /2 f 2 ).
See Table Supplementary Table 1 for the conversion between the desired group velocity ranges and the required value of the parameter B.

F. Synthesizing pulsed Bessel beams with separable spatio-temporal spectrum
To compare the performance of 3D ST wave packets with that of the conventional pulsed Bessel beams shown in Fig. 5(a-b) in the main text, we bypass the spectral analysis and 1D spectral transformation stages and send the input laser pulses directly to the 2D transformation.
Consequently, separable pulsed Bessel beams are produced. In addition, we spectrally filter ∆λ≈0.3 nm of the input pulses to obtain a spectral bandwidth comparable to that of the 3D ST wave packets.

Positive-luminal v=c B=∞
Positive-superluminal c< v<∞  The axial evolution of the time-averaged intensity profile I(x,y,z)= dt|ψ(x,y,z;t)| 2 is obtained in the physical space by scanning CCD 1 (The ImagingSource, DMK 27BUP031) along the propagation axis z [ Supplementary Fig. 8(b)]. This measurement confirms the diffraction-free propagation of the 3D ST wave packets, as shown in the Fig. 5 in the main text.

C. Time-resolved intensity measurements
The spatio-temporal intensity profile of the 3D ST wave packet is measured using a Mach-Zehnder interferometer in which the spatio-temporal synthesis arrangement is placed in one arm. The second (reference) arm is traversed by the original short plane-wave laser pulses from the source (pulse width ∼100 fs) that encounter an optical delay line τ [ Supplementary Fig. 16].
The two wave packets (the 3D ST wave packet and the reference pulse) propagate co-linearly after they are merged at the beam splitter BS 2 . When the two wave packets overlap in space and Supplementary Fig. 16. Schematic illustration of the setup for reconstructing the spatio-temporal profile of 3D ST wave packets and estimating their group velocity. The spatio-temporal synthesis arrangement from Supplementary Fig. 8(b) is placed in one arm of a Mach-Zehnder interferometer, and an optical delay τ is placed in the reference arm that is traversed by the short pulses from the initial laser source.
time, spatially resolved interference fringes are recorded by CCD 1 at a fixed axial position z, from whose visibility we extract the spatio-temporal profile at that plane. By scanning the delay τ, we can reconstruct the spatio-temporal intensity profile I(x=0,y,z;τ) at the fixed axial plane z.
Using this approach, we obtained the wave packet profiles shown in Fig. 6(b-d) in the main text (see [3,28] for further details).
To estimate the group velocity v of the 3D ST wave packet, we first arrange for the 3D ST wave packet and the reference pulse to overlap at z=0 as described above. We then axially displace CCD 1 to a different axial position ∆z, which results in a loss of interference visibility due to the group-velocity mismatch between the 3D ST wave packet traveling at v=ctanθ (where θ is the spectral tilt angle) and the reference pulse traveling at v=c. The interference is retrieved, however, by adjusting the delay by ∆t in the reference arm, from which we obtain the relative group delay between the two wave packets and thence the group velocity v= ∆z ∆t for the 3D ST wave packet [ Fig. 6(e) in the main text]. By repeating this procedure for ST wave packets with different θ (by varying the parameter B), we obtain the data plotted in Fig. 6(f) in the main text.
The uncertainty in the group-velocity measurement is estimated using the propagation-oferrors principle [29]. In our case, the largest contribution to errors in estimating v stems from the uncertainty δt in estimating the group delay ∆t, which is limited by the pulse-width ∆T, Supplementary Fig. 17. Schematic depiction of the off-axis holography strategy for reconstructing the complex field for 3D ST wave packets. (a) The system is similar to that in Supplementary Fig. 16 except that the reference pulse is spatially tilted with respect to the 3D ST wave packet. (b) Outline of the procedure for estimating the field amplitude and phase from the spatially resolved interference fringes.
which we set at δt=∆T/10. The error in the estimated value of v is δ v=| ∂ v ∂t |δt= v 2 ∆z δt. Using this relationship, we calculate the error bars in Fig. 6(f) in the main text.

D. Measurements of the field amplitude and phase for 3D ST wave packets
Finally, we make use of off-axis digital holography [30][31][32] to obtain the field amplitude and phase of 3D ST wave packets at fixed locations along z and at fixed instances in time τ. This is especially crucial to measure the phase of the OAM-carrying 3D ST wave packets shown in Fig. 7(b) in the main text.
We make use of the off-axis digital holography (ODH) methodology [30][31][32] to reconstruct the complex field of 3D ST wave packets ψ(x,y,z;t)=|ψ(x,y,z;t)|e iφ(x,y,z;t) . For this purpose, the same Mach-Zehnder configuration from the previous section is exploited but with a slight modification Supplementary Fig. 18. Measured phase profiles for 3D ST wave packets with and without OAM. (a) Measured transverse phase profile φ(x,y) (first row) and on-axis phase profile φ(x,y=0) along the x-axis (second row) at delays τ ≈0,±5 ps for a superluminal 3D ST wave packet. Here z=30 mm and =0; i.e., there is no OAM structure. (b) Same as (a) except that =1; i.e., OAM structure has been introduced into the field.
-we add a small angle between the propagation directions of the short reference pulse and the 3D ST wave packets, which are arranged to overlap in space and time at a fixed axial position z [ Supplementary Fig. 17(a)]. The interference pattern captured on CCD 1 contains a constant background term and an interference term of interest [ Supplementary Fig. 17(b)]. Following the ODH algorithm, we perform a digital fast Fourier transform (FFT) to separate the constant background from the interference term. By digitally isolating the first diffraction order of the Fourier-transformed image, we access the interference term that contains the complex field.
Finally, the inverse FFT of the centered first term gives the amplitude |ψ(x,y,z;τ)| and phase φ{ψ(x,y,z;τ)} of the 3D ST wave packets at that axial location z and time τ (see [30][31][32] for more details). By repeating this procedure for several instances τ we obtain the data plotted in Fig. 7 of the Main text.
The measured phase structure is plotted in Supplementary Fig. 18(a) for =0 and in Supplementary Fig. 18(b) for =. In the first case, in absence of OAM, the phase distribution takes a Gaussian form. The curvature increases as we move away from the center of the wave packet τ=0. We highlight the phase along the center of the wave front y=0, and especially where the intensity of the field is appreciable (between the two vertical dotted lines). Here the phase is flat at τ=0, and becomes Gaussian as τ increases. This is to be expected for a monochromatic Gaussian beam going through its focal point (where the phase front is flat). However, this is seen here in the time domain, which is a manifestation of so-called 'time diffraction' [33][34][35].
In presence of non-zero OAM, superimposed on the above-described behavior of the phase front is a helical phase distribution. At τ=0 the helical phase front flattens out, but appears as we move away from the wave packet center [ Supplementary Fig. 18(b)]. Once again, this is the expected behavior in space for a diffracting OAM mode when examined from the beam waist outward, but observed here in the time domain.